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ABSTRACT 

The Kustaanheimo-Stiefel transform turns a gravitational two-body problem into a 
harmonic oscillator, by going to four dimensions. In addition to the mathematical- 
physics interest, the KS transform has proved very useful in iV-body simulations, 
where it helps handle close encounters. Yet the formalism remains somewhat arcane, 
with the role of the extra dimension being especially mysterious. This paper shows 
how the basic transformation can be interpreted as a rotation in three dimensions. 
For example, if we slew a telescope from zenith to a chosen star in one rotation, 
we can think of the rotation axis and angle as the KS transform of the star. The 
non-uniqueness of the rotation axis encodes the extra dimension. This geometrical 
interpretation becomes evident on writing KS transforms in quaternion form, which 
also helps derive concise expressions for regularized equations of motion. 
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1 INTRODUCTION 

The Kustaanheimo-Stiefel transform is a remarkable rela- 
tion between the two most important elementary problems 
in dynamics: under a transformation of coordinates and 
time, a Kepler problem changes into a harmonic oscillator. 
Especially noteworthy is that the collision singularity in the 
Kepler problem is transformed into a regular point. The 
name comes from the works by Kustaanheimo (1964) and 
Kustaanheimo & Stiefel (1965), while the book by Stiefel & 
Scheifele (1971), which is largely devoted to the KS trans- 
form and its consequences, is perhaps the best known source. 
For a very short summary, see 'regularization' in Binney 
& Tremaine (2008). An important application of the KS 
transformation is in numerical orbit integration, where the 
singularity-removal is used to great advantage for simulating 
dense stellar systems with near collisions (Aarseth & Zare 
1974a,b; Mikkola & Aarseth 1990, 1993; Jernigan & Porter 
1989). Some recent papers also re-examine the formalism 
itself (Bartsch 2003; Waldvogel 2006). 

In two dimensions there is a much simpler version of 
the KS transform going back to Levi-Civita (1920). In the 
Levi-Civita transform, the coordinate plane is read as the 
complex plane, and the complex square root of the coordi- 
nate becomes the transformed coordinate. The geometrical 
interpretation is clear: the complex phase gets halved. The 
KS transform is also a kind of square root, but in four di- 
mensions. One wonders how the geometrical interpretation 
generalizes. 



It turns out that slewing a telescope is a convenient geo- 
metrical analogy. Suppose the telescope is at zenith, and we 
want to slew it to a particular star in one rotation. Normally 
we would simply move along a great circle from the zenith to 
the star. But we might prefer a different rotation (to avoid 
crossing the moon, say). For example, we could choose the 
midpoint on the above great circle and rotate about it by 
180°. In any case, the chosen rotation axis and the rotation 
angle are effectively the KS transform of the star. This idea 
of rotation in three dimensions about a non-unique axis gen- 
eralizes the idea of halving the phase in a complex square 
root. 

This paper attempts to provide some new insight into 
the KS transform by providing some reformulations and new 
derivations of known results, and especially to make the ge- 
ometrical interpretation evident. 



2 QUATERNIONS AND ROTATION 

Before considering KS theory, it is useful to review a concise 
algebraic way of specifying rotations in three dimensions, 
not often used in astrophysics but standard in computer 
graphics: quaternions. 

Quaternions are a generalization of complex numbers. 
The V— 1 °f complex numbers is replaced by three unit 
quaternions i, j, k, such that 

ii = jj — kk = —1, ijk = —1. (1) 
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From (1) it follows that ij = k = —ji and so on. In other 
words, quaternions are like a combination of dot and cross 
products in vector algebra. (Although historically quater- 
nions came before, having been invented by none other than 
W.R. Hamilton of Hamilton's equations.) 
A general quaternion has the form 



denote a point in space. The KS transform of q is the quater- 



A = A + A x i + A y j + A z k 



(2) 



where we will call Ao the real part. A quaternion with no 
real part is effectively a vector in three dimensions. 

In analogy with complex numbers, we will use the fol- 
lowing notation for quaternion conjugates and absolute val- 
ues. 



re [A] 
A* 



Ao 
A - 



It is easy to see that re [A*] = 
and as a result 

re [AB] = re [B* A*] = re [BA] 



Aii - A 2 j - A 3 k 
= Al + A\ + Al + A\ 

re [A] and (AB)* 



(3) 



B*A* 



(4) 



Rotation in quaternion notation is beautifully concise. 
Say we want to rotate a vector r by angle w about a unit 
vector n. Using quaternion algebra the rotation is simply 



R*rR 

where 

R — cos iw + sin \uj n. 



(5) 



(6) 



Unlike the equivalent expression using Euler angles, the 
expression (5) has no coordinate singularities (or "gimbal 
lock") and as a result is numerically more stable, which ex- 
plains its popularity in computer graphics. 

For an arbitrary (i.e., non-unit) quaternion R, the ex- 
pression (5) amounts to a rotation combined with scalar 
multiplication. 

It is possible to represent quaternions as matrices 
(though not necessary, even for numerical work). A familiar 
representation is in terms of Pauli matrices 



i = iai j — —i<72 k = iaz 



i 

1 



k 



(7) 



(8) 



Pauli matrices are most important as operators on quan- 
tum two-state systems (being Hermitian, whereas quater- 
nions are anti-Hermitian) . In recent years the most excit- 
ing two-state quantum systems have been Qbits in quantum 
computing. It turns out that expressions of the type (5) 
appear in the description of quantum-computing gates (see 
Mermin 2007, who also provides a derivation of essentially 
the above three-dimensional rotation formula). 



3 THE KUSTAANHEIMO-STIEFEL 
TRANSFORM 



Let 

q = xi + yj + zk 



Q = Qo + Qxi + Q y j + Qzk 
the transformation formula being 
q = Q*kQ. 
A solution for Q is 
_ xi + yj + Zk 



>2Z 



Z = z+^x 2 +y 2 + z 2 



(10) 



(11) 



(12) 



as is easily verified by multiplication, following the quater- 
nion rules. But Q 1 is not unique, because changing to 



Q = (cos ip — sin ip k) Q 1 



(13) 



leaves Equation (11) invariant. Thus ip behaves like a gauge. 

Everything so far is already in the literature. The new 
result in this paper is that we can readily visualize Q, in- 
cluding its non-uniqueness. 

Comparing (11) and (5), it is evident that Q is a rotator 
that takes the z axis to q. To visualize Q, let us rewrite q 
as 



q = r(sin 9 cos (pi + sin 6 sintpj + cos 8 k) 



(14) 



where r, 9, <f> are the usual polar coordinates. Rewriting Q 1 
in the solution (12) and simplifying, we have 

Q 1 = v^sin ^6 cos (pi + sin ^0sm</>j + cos \0k). (15) 

In other words, the zenith distance of Q 1 is halfway along the 
great circle from k to q. From (6) we see the rotation angle 
lo would be 7r. Now let us apply the gauge transformation 
(13) with ip = it/ 2 to Q 1 . This gives 

Q 11 = v^cos |# + sin ^9 sin (j>i — sin |0 cos 4>j) (16) 

Now the implied rotation is by 9, about an axis perpendic- 
ular to both k and q. In general, we can write 



Q = cos ip Q 1 — sin ip Q 1 



(17) 



(9) 



which is to say, Q could be anywhere on the great circle 
joining Q 1 and Q 11 . The telescope-slewing analogy given 
above is simply a description of the preceding three formulas. 

An interesting special case is <p> = 0, which gives q — 
r(cos9k + sinOi) and Q 1 — ^/r(cos \9k + sin \9i). Then 
Q 1 is effectively the complex square root of q (we need to 
read k as the real axis and i as the imaginary axis). In other 
words, the planar case can be reduced to the Levi-Civita 
transform by a suitable gauge. 

Quaternion formulations of the KS transform have been 
discussed in several authors: Stiefel & Scheifele (1971) men- 
tion quaternions but appear to dislike them, while later au- 
thors (for example Vivarelli 1994; Waldvogel 2006) are more 
favourable. The precise definition adopted for the transform 
varies, but is equivalent to Eq. (11). That Q represents a ro- 
tation and shrinking/stretching of q is also known. Bartsch 
(2003) specifically notes that the rotation axis is unique in 
two dimensions but not in three. But the explicit description 
of the implied rotations, as above, appears to be new. 
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4 THE CANONICAL MOMENTUM 

So far we have just discussed geometry, but of course the 
real significance of the KS transform is dynamics, which we 
now consider. Let 

p = Pxi + Pyj + Pzk (18) 

be the canonical momentum conjugate to q. We seek 

P = Po + Pxi + Pyj + P*k (19) 

that will be canonically conjugate to Q. Let us write 

re [p*dq] = re \p* dQ* kQ] + re [p*Q*k dQ] (20) 

Using the identity (4) we can rewrite the middle term 
as re \p (dQ* k Q)*]. Since p — —p* and (dQ* kQ)* — 
Q* (— k) dQ the term becomes becomes re [p*Q*kdQ]. Thus 
we have 



ro[p*dq] — 2rc [p*Q*kdQ] 
Now if we define 
P = -2kQp 

we have 

rc [p*dq] — re [P* dQ] 



(21) 



(22) 



(23) 



which is to say, p ■ dq = P ■ dQ. Provided the Hamiltonian 
depends on P, Q only through p, q and not on the gauge ip, 
the transformation (P,Q) — > (p,q) is canonical. 

To get an explicit expression for p, we multiply (22) on 
the left by Q*k, obtaining 



P 



Q*kP 
2Q 2 ' 



(24) 



Note that while we have to be careful about the order of 
multiplication when i,j,k are involved, real numbers like 
Q 2 commute with everything. Since p has no real part, 
rc [Q*kP] = identically. We can think of it as a formal 
constant of motion resulting from invariance with respect to 
i>. 

That P (as defined in Equation 22, or equivalently) 
completes a canonical transformation is a standard part of 
KS theory, but the derivation of the canonical condition us- 
ing quaternion identities appears to be new. 



5 THE TWO-BODY PROBLEM AND THE 
HARMONIC OSCILLATOR 

Let us now write the Kepler Hamiltonian 



H = ip 2 - 1/9 



(25) 



in terms of KS variables. Multiplying each of (11) and (22) 
by its quaternion conjugate, we have 



2 ^.4 

q =Q , 



Ap 2 Q 2 



and substituting these gives 
H =\P 2 /Q 2 -1/Q 2 . 



(26) 



(27) 



We now use a device known in Hamiltonian dynamics as a 
Poincare time transformation. This involves introducing a 
fictitious time variable s, whose relation to t we choose to 
be 



dt = Q 2 ds. (28) 

Since Q 2 is the radial distance in the Kepler problem, (28) 
is in fact Kepler's equation, and s is the eccentric anomaly. 
In the fictitious time variable s, the equations of motion are 
given by a new Hamiltonian 



T = Q (H - E) 



I p 2 

8 



EQ 2 



1 



(29) 



with E being the constant initial value of H. The time- 
transformed T Hamiltonian is zero along a trajectory, but 
its partial derivatives are not zero. 

The Hamiltonian V is remarkable indeed. For E < 
(bound orbits) it is a harmonic oscillator. Since Q has four 
components, F is like a mass on an isotropic spring in four 
Euclidean dimensions. Thus the well-known fact that the 
bound Kepler problem has a dynamical 0(4) symmetry. 
For the unbound case, the symmetry group is different: for- 
mally the Lorentz group, but with a physical meaning com- 
pletely different from special relativity. And — perhaps most 
importantly — Hamilton's equations for T are well-behaved 
even at Q — (a collision). This is known as regularization 
and was the original motivation for KS theory. 

The effect of an external force F is simple. From (22) it 
follows immediately that F will add an extra contribution 
of —2kQF to dP /dt, which amounts to a contribution of 
—2Q 2 kQF to dP/ds. Provided the external force is non- 
singular, the equations of motion in s remain regular. 



6 REGULARIZING THE THREE-BODY 
PROBLEM 

Application of KS regularization to iV-body simulations in- 
volve expressing the gravitating system either as a tree-like 
hierarchy of coupled two-body systems (Jernigan & Porter 
1989) or as a chain (Mikkola & Aarseth 1990, 1993). The ba- 
sic idea can be described using the three-body problem with 
all masses unity. Here again, quaternions enable a concise 
formulation. 

In relative coordinates, the Hamiltonian for three unit 
gravitating masses can be written (cf. Eq. 12 in Aarseth & 
Zare 1974a) as 

H=\p\ + \pl+ Pl -p 2 ---- (30) 

qi 92 

plus an additional potential V(q 1: q 2 ). Here q x ,q 2 expresses 
the position of the first and second body relative to the ze- 
roth body, while p lt p 2 express the momenta of the first 
and second bodies in the barycentric frame. Meanwhile, 
V(q 1 ,q 2 ) expresses the mutual interaction of the first and 
second bodies, plus any external potential. We can regard 
V(q 1 ,q 2 ) as an external potential, and since we already 
know how to deal with external forces, we set V aside and 
concentrate on H. 

Now we introduce KS variables q 1 = Q{kQ 1 and so on. 
Defining 



n = rc [(Q* 1 kP 1 )*Q* 2 kP 2 



(31) 



we can write p 1 ■ p 2 as H/(4Q 2 Q 2 ). Applying a Poincare 
time transformation 



dt = Q\Q%ds 
gives 



F = Q\Q\(H - E) 



(32) 
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r = |Pi 2 Q? + \PlQl + |n - Ql - Qi - EQlQl (33) 

where E is the value of H. The F Hamiltonian has no de- 
nominators, and is thus regular for collisions with the zeroth 
body. (We assume V remains regular, that is to say, the first 
and second bodies do not collide with each other or any other 
bodies than the zeroth. In practice, simulations redefine the 
relative coordinates whenever necessary, according to who 
is close to whom.) 

For the equations of motion we need derivatives with 
respect to quaternion components. First we have Vq 1 Q^ = 
2Q 1 . Slightly more subtle is Vq 1 re [Ql A] = A if A is in- 
dependent of Q 1 . Using this last identity, together with the 
definition (31) of n, we derive 

V Pi n = kQ L Q* 2 kP 2 , V Ql II = kP 1 P* 2 kQ 2 . (34) 
Hamilton's equations are then 

^gi = 1 1 Q(p 1 -kP 1 p;kQ 2 

= (2 + 2EQl-\P*)Q 1 +kQ 1 QlkP 2 (35) 
and similarly for P 2 ,Q 2 . 

7 DISCUSSION 

In dynamical astronomy the KS transformation is profound, 
but may appear mysterious. This paper attempts to make it 
less mysterious, and hopefully therefore more useful, by ex- 
plaining it in three-dimensional geometric terms. There are 
several possible directions in which the KS transformation 
may turn out to be useful. 

First, one can imagine new orbit integrators specialized 
to nearly-Keplerian problems. Work on dense stellar sys- 
tems with near collisions has already been mentioned (for 
reviews see the books Aarseth 2003; Heggie & Hut 2003). In 
the planetary regime, which differs from the dense-stellar 
case in having few bodies but many more orbital times, 
time transformations reminiscent of (28) used for KS reg- 
ularization have proved useful for highly eccentric orbits 
(Mikkola 1997; Emel'yancnko 2002), while some integration 
algorithms (Mikkola & Tanikawa 1999; Preto & Tremainc 

1999) apply the time transformation (28) implicitly. Could 
the KS transformation itself be exploited here? Fukushima 
(2005) has some further ideas. 

Second, it is conceivable that KS variables could sim- 
plify perturbation theory. Perturbation theory in classical 
celestial mechanics (see for example Murray & Dermott 

2000) is algebraically frighteningly complicated, basically 
because the natural variables for the unperturbed and per- 
turbed parts (being the Keplerian action-angles and real- 
space coordinate) are related through an implicit equation. 
On the other hand, the action-angles of the KS-transformed 
Kepler problem are explicitly related to space coordinates — 
the implicit equation is transferred to the time variable. 
Could some major simplication be achieved through KS vari- 
ables? Some progress has been made by Vrbik (2006). 

Third, the KS transformation might provide new in- 
sight into analogous quantum problem. Bander & Itzykson 
(1966a,b) derive the symmetry groups of the bound and un- 
bound Coulomb problem. These turn out to be the same 



four-dimensional symmetries as in KS theory. Is the KS 
transformation implicit in that work? 
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